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Abstract 

The response of the solar coronal magnetic field to small-scale photospheric boundary motions in- 
cluding the possible formation of current sheets via the Parker scenario is one of open questions of solar 
physics. Here we address the problem via a numerical simulation. The three-dimensional evolution of 
a braided magnetic field which is initially close to a force-free state is followed using a resistive MHD 
code. A long-wavelength instability takes place and leads to the formation of two thin current layers. 
Magnetic reconnection occurs across the current sheets with three-dimensional features shown, including 
an elliptic magnetic field structure about the reconnection site, and results in an untwisting of the global 
field structure. 
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Figure 1: TRACE coronal loops. (Left) A large-scale tangled configuration and {right) apparently smooth 
loops. 

1 Introduction 

Parker's notion of the 'topological dissipation' of coronal magnetic fields (Parker, 1972) continues to 
generate much debate. Simply put, Parker's suggestion is that following boundary motions of sufficient 
complexity the magnetic field of a coronal loop will be unable to ideally relax to a smooth force-free equi- 
librium and instead tangential discontinuities in the field, corresponding to current sheets, will develop. 
In general terms the possible outcomes of relaxation are a development of singular current sheets (e.g. 
Ng & Bhattacharjee, 1998; Janse & Low, 2009), development of thin but non-singular current layers (e.g. 
Longcope & Strauss, 1994; Galsgaard et al. 2003), and a smooth equilibrium with large-scale current 
features (e.g. van Ballegooijen 1985; Craig & Sneyd, 2005, Wilmot-Smith et al. 2009a). It should also 
be noted that the distinction between the first two cases may be difficult to determine numerically. 

In a previous a paper (Wilmot-Smith et al. 2009a) we considered the ideal relaxation of a braided 
magnetic field towards a force- free equilibrium. The magnetic field configuration was based on the pigtail 
braid and imposed as an initial condition (rather than being built up through boundary motions) . This 
particular braid was chosen since it is the simplest non-trivial braid with no net twist and, accordingly, 
is the most conservative realistic case to examine. More complex braided fields will, in general, con- 
tain components of the pigtail-type. There is ample motivation for modelling loops as having braided 
components. Photospheric turbulence subjects loop footpoints to a random walk, with motions of the 
fragments about each other acting to braid (or unbraid) the overlying loop. However, while there is 
some evidence for the existence of braided loop configurations (see the left-hand image of Figure [l]for 
example), most coronal loops appear to be close to a potential field (e.g. Figure [l] right-hand image). 
This may be an effect of the large aspect ratios typical to loops, since a winding of a field line around 
another field line is almost undetectable when smoothed out over the length of the loop (see also Berger 
& Asgari-Targhi, 2009). Another reason could be that reconnection is very efficient in maintaining a low 
degree of topological complexity in loops. The present work is designed in part to test whether this is 
the case. 
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In Wilmot- Smith et al. (2009a) an ideal Lagrangian relaxation scheme (Craig & Sneyd 1986; Pontin 
et al. 2009) was used to ideally evolve the braided field towards a force-free state. A smooth near force- 
free equilibrium was attained with large-scale current features i.e. without any tangential discontinuities 
or even strong current concentrations. (This situation may be compared with Parker (1994) where the 
same pigtail braid situation is considered as a thought experiment and concluded to inevitably lead to 
tangential discontinuities. It appears that Parker's assumption of a = in the domain fails; we indeed 
have / a (iS* = 0, where 5* is a cross-sectional surface through the braid, but a varies significantly between 
field lines). 

Although the local current density in the ideal equilibrium is of large-scale, a global quantity, the 
integrated parallel current, J J\\ dl, has small scales (here the parallel indicates parallel to the magnetic 
field and the integral is taken along magnetic field lines). By small scales we mean that f J\\ dl varies 
significantly between neighbouring field lines as the field lines pass though a particular plane, such as the 
lower boundary. In Wilmot-Smith et al. 2009a it was suggested that these small scales in J J\\ dl could 
lead to the development of a resistive instability and so a loss of equilibrium of the field. This mechanism 
would be distinct to that of Parker's topological dissipation but have the same consequence. 

The consideration that for sufficiently small scales in the integrated parallel current then a force- 
free magnetic field will be resistively unstable (the instability being dependent, for a given scale in the 
integrated parallel current, on the value of the resistivity) motivated us to put the end, near force-free, 
state of the Lagrangian relaxation into a resistive MHD code and test whether or not the state is stable 
(for various values of resistivity allowed by numerical limitations). It turns out that the braided field as 
implemented in the resistive code is not stable. This finding, together with the subsequent evolution of 
the field, allows us to address a number of key questions related to MHD and the behaviour of coronal 
magnetic fields. The results are described in this series of papers. Here we describe the details relating to 
the numerical setup of the problem and the early evolution of the system. A subsequent paper addresses 
the long term evolution of the system. 



2 Numerical Scheme and Simulation Setup 
2.1 Numerical Scheme 

The numerical scheme employed in the simulations that follow is described briefiy below (further details 



may be found in Nordlund & Galsgaard (1997) and at http://www.astro.ku.dk/~kg). We solve the 
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three-dimensional resistive MHD equations in the form 



E = - (v X B) + ryJ, (2) 

J = V X B, (3) 

I = -V.(pv), (4) 

^ (pv) = -V • (pvv + r) - VP + J X B, (5) 

1^ = -V • (ev) - P V • V + Q^i,, + Qj, (6) 

where B is the magnetic field, E the electric field, v the plasma velocity, rj the resistivity, J the electric 
current, p the density, r the viscous stress tensor, P the pressure, e the internal energy, Qvisc the viscous 
dissipation and Qj the Joule dissipation. An ideal gas is assumed, and hence P — (7 — 1) e = |e. 
These equations have been non-dimensionalised by setting the magnetic permeability /io = 1, and the 
gas constant (7^) equal to the mean molecular weight (M). The result is that for a volume in which 
\p\ — |B| = 1, time units are such that an Alfven wave would travel one space unit in one unit of time. 

The equations (1-6) are solved on staggered meshes; with respect to a mesh on which p and e are 
defined at the body centre of the cell, B and P are defined at face centres and E and J at edge centres. 
In this way the required MHD conservation laws are automatically satisfied. Derivatives are calculated 
using sixth-order finite differences that return a value which is shifted half a grid-point up or down with 
respect to the input values. When the staggered mesh means that some quantity must be interpolated, 
data values are calculated using a fifth-order interpolation method at the relevant position. A third-order 
predictor-corrector method is employed for time-stepping. 

In all simulation runs we employ a spatially uniform resistivity model. Viscosity is calculated using a 
combined second-order and fourth-order method (sometimes termed 'hyper- viscosity'), which is capable 
of providing sufficient localised dissipation where necessary to handle the development of numerical 
instabilities (Nordlund & Galsgaard 1997). The effect is to 'switch on' the viscosity where very short 
length scales develop, while maintaining a minimal amount of viscous dissipation where the velocity field 
is smooth. 



2.2 Creating the Initial Condition 

As discussed in Sec. ^ the initial state for the resistive MHD simulations is drawn from the final state 
of the Lagrangian relaxation experiment of Wilmot-Smith et al. 2009a. The quantities previously known 
from the Lagrangian code (see Craig & Sneyd 1986) are the magnetic field B and the current J and in 
the near force-free relaxed state these are known on a highly distorted mesh. We describe below the 
process of constructing the field on the rectangular grid required for the resistive MHD simulations. 
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In order to ensure that the initial magnetic field is divergence- free, we work first with the vector 
potential A for B. In the relaxation scheme, the calculation of A requires only that we know the initial 
vector potential before relaxation and the mesh deformation. In terms of the initial mesh X and the final 
'relaxed mesh' x, the ith component of the final vector potential A-^ is given (see Appendix) in terms of 
the initial vector potential A° by 

To create the input magnetic field for the MHD simulations we then interpolate this vector potential 
onto a rectangular grid. Since the magnetic field components are face-centred on the staggered grid, the 
vector potential components are interpolated onto locations corresponding to edge centres of the desired 
grid. We then obtain the magnetic field by taking the curl of A using the sixth-order finite differences 
described above, which yields magnetic field components at face-centres as required. An interpolation 
scheme using biharmonic spline radial basis functions was applied to A, the particular scheme chosen to 
maximize the smoothness of the corresponding current density J, which involves second derivatives of 
A. To further improve this smoothness a simple five-point smoothing algorithm was finally applied to 
A, before taking the curl. 

The result of the above is that the initial braided magnetic field for our MHD simulations is divergence- 
free to accuracies on the order of truncation errors of the sixth-order finite differences (with typical 
maximum |V • B| ^ 10~^ within the domain). The topology of the magnetic field turns out to be well 
conserved by the process, another important consideration for the experiment. However, a drawback is 
that the quality of the force-free approximation is not perfectly maintained; the initial state is further 
from force-free than the relaxed field of the Lagrangian experiment. Details and implications are discussed 
in Sees. [3] and m 

2.3 Initial State 

The initial magnetic field is given on a domain of size [—24, 24] in the vertical (z) direction and [—6, 6] 
in both the horizontal (x, y) directions. Covering this domain we take a uniform mesh of 320"^ cells 
and employ closed boundary conditions in all three directions. The magnetic field is line-tied, and can 
be very closely approximated by B = [0, 0, 1] at the boundaries, i.e. it is directed perpendicular to 
the z boundaries and parallel to the x and y boundaries. To achieve the perpendicular condition the 
Lagrangian relaxation experiment was re-run on the larger horizontal domain (now [—6,6]^ compared 
with [—4,4]^ in Wilmot-Smith et al 2009a). The braided field is centered in the middle of the domain 
and the field is close to uniform in the region external to the braid. Accordingly the results presented here 
are shown for only the subsection of the full domain in which the important dynamics occur, specifically 
[-3,3]2 X [-24,24]. 

An isosurface of current in the initial state is given the left-hand image of Fig. [2] The current has 
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Figure 2: Initial state of the simulation: {left) Isosurface of current |J| at 25% of the domain maximum and 
(right) some particular magnetic field lines illustrating the braided nature of the field. 



large-scales (see also the upper-left hand image of Fig. |5] where contours of current in a horizontal cross- 
section are shown) with two fingers of current extending vertically through the domain. Some sample 
field lines further illustrating the nature of the field are given in the right-hand image of the same figure. 
Although non-trivial, the initial state has little magnetic energy in excess of the homogeneous field (0.96% 
in [—3, 3]^ X [—24, 24]). The aspect ratio employed in the model is 1:8. Although this is larger than that 
of many previous simulations it is smaller than a realistic ratio for a coronal loop (1:50, say). In the 
configuration the poloidal field components are small compared with the toroidal components so that the 
field lines look almost straight. This level of braiding would be observationally difficult to distinguish 
from a potential field. 



To initialise the simulation the dimensionless plasma density (Sec. 2.1) is set at p = 1 throughout 
the domain and the internal energy as e = 0.1. The result is a plasma- /3 at t = that lies in the range 
/3 G [0.1,0.14]. For the results described in the main section of this paper (Sec.|3| we consider the early 
evolution of the system (up to t = 14) with time measured in units of the Alfven time. A uniform 
resistivity of 77 = 0.001 has been taken and the effect of changing the resistivity is discussed at various 
points in the following text. 
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3 Results 



3.1 Basic Properties 





Figure 3: Maximum absolute value of the current {top left), velocity {top right) and Lorentz force {bottom 
left) and total kinetic energy {bottom right) in the domain with time for a sequence of decreasing uniform 
resistivies as indicated in the figures. 

Figure [S] shows the maximum absolute values of the current, velocity and Lorentz force, and the 
kinetic energy (J ^v^dV) for the time interval (t G [0, 14]) under consideration. The domain taken in all 
cases is the central section, [—3,3]^ x [—24,24], of the full box. The variation in quantities is shown for 
a sequence of uniform resistivities decreasing by over an order of magnitude, specifically 77 = 0.005, 0.001 
and T] = 0.0002. 

The initially high maximum Lorentz force, |j x B|macc = 0.501 decreases rapidly over the first few 
time units. Both the high value and the decrease are artifacts of the method used to create the initial 



state. The interpolation required to transfer the state to the Eulerian grid (as described in Sec. 2.2) 
results in some noise in the initial magnetic field and current density. Some noise persists even with the 
application of a smoothing algorithm to the vector potential for the magnetic field and this is particularly 
noticeable in the Lorentz force rather than the magnetic field and current density alone. Thus whilst the 
final state of the Lagrangian relaxation experiment had a very small maximum Lorentz force (specifically 
|j X B I mace ~ 2 X 10~^), the initial state here is further from force- free. The decrease over t G [0, 2] then 
arises though a smoothing of the noise in the initial state. 

Turning now to the remaining quantities shown in Fig. [3] two primary features are found. Firstly a 
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growth in the kinetic energy and maximum velocity for t G [0, 4] occurs, and this growth is independent 
of resistivity, 77. For the remaining time considered there is no significant change in kinetic energy. 
Secondly, a linear growth in the current density occurs for t G [8, 12]. The rate at which the maximum 
current density increases is higher for lower resistivity, 77. At t = 12 the maximum current density is 
achieved; this maximum is higher for lower resistivity but for all three resistivities the growth phase ends 
at the same time. 

The lack of dependence of kinetic energy on resistivity may suggest an ideal instability is present. 
The subsequent linear growth of current would then be a non-linear consequence of this instability rather 
than its initial appearance. This growth is clearly dependent on resistivity, being slower for higher values 
of 77. Little is known about the no n- linear phase of instabilities and such a dependence may still be 
consistent with an ideal instability with a no n- linear phase damped by ry. Strong conclusions are clearly 
difficult to draw at this stage. An additional consideration is that the implementation of the field on 
the new grid has resulted in significant Lorentz forces in the initial state. We return to these questions 
in Section [4] but now proceed to consider the the nature of the currents within the domain, now fixing 
r] = 0.001. 

3.2 Formation of Current Layers 

Figure [4] shows isosurfaces of current at 50% of the domain maximum ( J — J max 

/2) for a sequence 

of increasing times (for the initial state see Fig. [2]). In the early stages (t G [0,4]) the current diffuses 
slightly while maintaining its large scales in the perpendicular direction. A symmetric evolution follows 
and after the phase of current growth two current concentrations are present, centered at 2; = 3.4 and 
z = —3.6. We call these the two 'initial current layers'. 

The formation of these two initial current layers is best illustrated by considering a horizontal cross 
section through the central plane = 0). Figure [s] shows contours of the vertical component of current 
{\Jz\) in that plane at t = 0,6, 12. The z-component is taken since it significantly dominates over the 
two horizontal components, as evident in the shape of the current sheets (Fig. [4]). Note that in order to 
incorporate both sheets the cross-sectional plane chosen, ^ = 0, does not intersect the centre of either 
current sheet and so the magnitude of current in this plane is somewhat low in comparison to the domain 
maximum. The collapse of the two oppositely signed large-scale fingers of current present in the initial 
state into two thin current sheets of correspondingly the same sign is clearly shown. Also evident is 
the formation of a weaker current envelope around the braided flux, separating it from the uniform 
background field. Cross-sections of \Jz\ in the horizontal planes through the centres of the two current 
sheets are shown in the final two images of Fig. |7| 
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Figure 4: Isosurfaces of current, |J|, at 50% of the domain maximum for a sequence of increasing times (from 
left to right, t = 4, 6, 8, 10, 12) showing the formation of the two initial current layers. 



t=0.00 t=6.00 




Figure 5: Contours of the vertical component of the current, J^, in the central plane, 2: = 0, at increasing 
times, illustrating the formation of the two initial current layers (first three images). The lower right-hand 
image shows contours of integrated parallel current along field lines in the initial state at the central plane 
z = 0. From this quantity field line tracers for the locations of initially high integrated parallel current have 
been determined and marked in crosses in the previous frames, as described in the main text. 



3.3 Predictors of Current Layers 

The field lines along which these two current sheets form turns out to be well predicted by the regions 
of high integrated parallel current in the initial state. In resistive MHD the integrated parallel current is 
related to the integrated parallel electric field via the relation / J^^dl = r] f E\\dl (in the case of a uniform 
resistivity 77, as considered here). The integrated parallel electric field is a key quantity for 3D magnetic 
reconnection; for a localised reconnection region the maximum value of / E\\dl determines the rate of 
reconnection. 

Shown in the lower right-hand panel of Fig. |5] are contours of the integrated parallel current, / J\\dl^ 
in the initial state, t = 0, where the path of the integral is taken over magnetic field lines. Again the 
contour is shown in the central plane (z = 0). To obtain this contour map, field lines have been integrated 
through 160^ grid of points covering the domain x,y ^ [—3,3]^, z = 0. Two peaks in the quantity are 
present and the structure is quite different to that of the current itself in the initial state. We now 
identify those field lines in this tracing procedure for which the value of | / J\\dl\ is greater than or equal 
to 75% of the domain maximum, noting the locations where they intersect the lower boundary (where 
the locations of the field lines are held fixed). For the sequence of times of Fig. |5] we trace field lines 
starting from those locations on the lower boundary up through the domain and mark with a cross in 
that same figure their point of intersection with the z — plane. Here the difference in colours indicates 
field lines with positive (black crosses) and negative (white crosses) integrated parallel current, although 
this distinction is made only to facilitate identification of the locations. It is found that these field lines, 
traced from the initial locations of high integrated parallel current, are good indicators for the locations 
of formation of the two current layers. Since the flux on the lower boundary is held fixed these may be 
identified as the same field lines for as long as the evolution remains ideal. Whilst the evolution will be 
ideal only during the early stages of the simulation it is clear that the tracers do, nevertheless, provide 
a useful predictor for the locations of current sheet formation. 

Locations of high integrated parallel current are not a commonly used indicator for current sheet 
formation. Indeed it is quasi-separatrix layers (QSLs), regions where the field-line connectivity varies 
strongly (Priest & Demoulin, 1995) that are widely thought of as likely sites of current sheet formation. To 
identify QSLs (as well as their intersections, hyperbolic flux tubes or HFTs) the squashing factor (Titov 
et al. 2002) is used. Usually denoted by Q, the squashing factor is an indictor of field line connectivity 
and takes on high values in regions where the field line mapping is strongly distorted. Regions of high 
Q outline QSLs. As discussed in Wilmot-Smith et al. (2009b), the braided magnetic field taken as the 
initial state here contains several QSLs. Contours of the squashing factor, Q, in the central plane {z — 0) 
are shown in Fig. [g] at t = (left) and t = 12 (right). For the calculation again 160^ points over the 
region x,y ^ [—3,3]^ have been used, a number comparable to the grid resolution. 

Several regions of high Q are present in both snapshots. The two regions of highest Q in the initial 
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t=12.00 




Figure 6: Contours of the squashing factor, specifically log^Q (Q), in the central plane, z = 0, for (left) the 
initial state, t = and (right) the time t = 12. 

state are associated with the two initial current layers at the later time, that is, the current sheets have 
formed along QSLs of the field. At the same time, several other regions of high Q are present (both at 
t = and t = 12) that are not associated with any particular current features. For example, in the initial 
state there are eight distinct regions where log^o (Q) is greater than 85% of its maximum value but only 
two of these regions correspond to particular current features at t = 12. These results suggest that the 
integrated parallel electric current and the squashing factor could be used in conjunction for predicting 
current sheet formation more accurately than by using Q alone. 



3.4 Plasma Flows and Reconnection 

We move now to consider the nature of the plasma flows within the domain. First recall that the low 
values of the plasma beta (Sec. [5]) imply the dynamics will be dominated by the Lorentz force rather than 
the gas pressure. With both the magnetic field and the current having stronger vertical than horizontal 
components we have that the Lorentz force is primarily in the horizontal direction and so, similarly, are 
the plasma flows. The Lorentz force drives a flow with a dipolar structure in the six regions of initially 
strongest current with direction dependent on the sign of the twist in that region. An illustration of such 
a flow is shown for the plane z = 3.4 (a negative twist region corresponding to one of the sites of current 
sheet formation) in the first (upper left) image of Fig. |7| 

The next four images in that figure show the development of the flow for a sequence of increasing 
times up to the point of maximum current. In the sequence the length of the arrows indicating flow 
direction have been normalized to each image. The background contours show the vertical component of 
current in that plane with the colour scale normalised to the current at t = 12 (lower-left- hand image), 
given by the first colour-bar. The sequence clearly shows the association of the stagnation part of the 
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dipolar flow with the current intensification. The out-flowing plasma from the location of stagnation sets 
up a counter-flow to the initial dipolar structure leading to oppositely directed flows on either side of the 
weak enclosing current sheath. The flnal image shows plasma flows and current in the plane z — —3.6 at 
t —VI. This cross- section is across the second current sheet and shows the naturally expected inversion 
of the flow direction. The result is that the global flow structure is dominated by rotational components 
the direction of which varies both vertically along the structure and on either side > 0, 2/ < 0) of the 
braided fleld. 

As already indicated by Fig.js] at t = 12 the maximum magnitude of the plasma flows (|V| max 0.24) 
is a signiflcant fraction of the Alfven speed (Va 1.2). These strongest flows are associated with the 
global rotation and not with outflows from the two initial current layers (which have an associated 
|V|macc ^ 0.15). However, we do expect magnetic reconnect ion to be taking place across these current 
sheets and so proceed to consider the nature of this reconnection. For this we concentrate again on 
the initial current layer centered at 2; = 3.4 (a similar situation occurs about the other current layer) 
and focus on the structure of the fleld and flows in the plane perpendicular to the magnetic fleld at the 
location of maximum current magnitude, | J|. 

In the left-hand image of Fig. [S] we indicate in more detail the nature of the stagnation flow about the 
current sheet. Superimposed on the background are contours of the out-of plane component of vorticity, 
i.e. the component of vorticity in the direction of the magnetic fleld at the location of maximum current. 
The vorticity shows a quadrupolar conflguration around the current sheet. 

Perhaps more of a surprise is the structure of the magnetic fleld in the region. The right-hand image 
of Fig. [8] shows the components of the magnetic fleld in the cross-sectional plane under consideration. 
The magnetic fleld is shown to have an elliptic conflguration about the current sheet, i.e. about the 
reconnection region. This flnding contrasts with the two-dimensional picture of magnetic reconnection 
under which the process can only occur at a hyperbolic (X-type) null-point of the magnetic fleld. In three- 
dimensions a much wider variety of possible reconnection sites exist. Reconnection may be associated 
with 3D null-points (e.g. Lau & Finn, 1990; Priest & Pontin, 2009), magnetic separators (which connect 
two null-points, e.g. Longcope & Cowley, 1996; Pontin & Craig, 2006; Haynes et aL, 2007), or may occur 
in the absence of any such topological features (Schindler et aL, 1988), the latter sometimes termed 'non- 
null reconnection'. In particular, the local magnetic fleld structure need not be hyperbolic but may be 
elliptic (Hornig & Priest 2003), as recently found in some 3D numerical simulations of reconnection. For 
example, Wilmot-Smith & De Moortel (2007) considered reconnection occurring along a quasi-separator 
and found an elliptic fleld structure in perpendicular cross-sections. The separator conflguration of Parnell 
et al. (2010) showed an elliptic structure along a signiflcant length of the separator under consideration. 
Parnell et al. (2010) also discussed the separator case theoretically, concluding an elliptic conflguration 
would be a generic situation given a suflficiently strong current density along the separator. Our flndings 
demonstrate that an elliptic field configuration may be present about reconnection sites in the non-null 
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t=2.00 



t=6.00 




ure 7: Arrows indicate plasma velocities [14, ^y] in horizontal cross-sectional planes, superimposed on 
vertical component of current, J^. The first five images show structures about the upper current sheet 

Lug the plane z = 3.4) while the final image is for the lower current sheet (taking the plane z = —3.6). 

images are given at various times, as indicated in the figures. 



13 




\'fy///. — 

^ \ \ \ \ V 
\ \ \ \ \ \ 
I I i 1 1 i 

^ / / / J / 
^ ^ ^ / / / 
^ ^ ^ ^ ^ / 
^ ^ ^^^^^^ ^ ^ 

^ ^ ^ 





Figure 8: Here various quantities are considered in the plane perpendicular to the magnetic field at the 
region of maximum (negative) current. {Left) Arrows indicate plasma flows in that plane while the coloured 
contours show the out-of-plane component of vorticity. (Right) Arrows indicate magnetic field components 
in the plane superimposed on contours of the out-of-plane current. The colour table for each is blue-yellow 
for negative - positive. 



case. Additionally, tracking these field lines back to the initial setup, an elliptic perpendicular field 
configuration is again found indicating that locally hyperbolic structures are not necessary for current 
intensification. As previously discussed, the squashing factor (Q) at the two reconnection sites is very 
high which demonstrates a further point; regions of highest-Q within a domain may have a locally elliptic 
field configuration. 



4 Nature of Instability 

Evidently the initial magnetic field configuration is not in a stable equilibrium. Since an exact equilibrium 
of the ideal relaxation code employed in Wilmot-Smith et ai (2009a) is known to be linearly ideally stable, 
the lack of stability could arise from one of a number of factors: 

1. A resistive instability. The relaxed state of Wilmot-Smith et ai (2009a) contains small scales in 
the integrated parallel current. Small enough scales (for a given resistivity) in this quantity are 
incompatible with a resistive equilibrium. 

2. A non-linear ideal instability not previously found in the ideal evolution of Wilmot-Smith et 
al. (2009a) since that evolution only guaranteed linear stability. 

3. A lack of equilibrium in the initial state either since: 

• The path to equilibrium of the ideal relaxation is fictitious and an exact equilibrium had not 
been reached. Numerical difficulties (described in Pontin et al. 2009) result in the final state 
of the ideal relaxation having |J x B|max ~ 2 x 10~^, i.e. is not perfectly force-free. Thus the 
final state of Wilmot-Smith et al. 2009a need not be stable (or, indeed accessible via a real 
MHD relaxation dynamics). 
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t=20, middle-third 




Figure 9: (Left) Vertical component of current in the central plane = 0) at t = 20 for the an MHD evolution 
on only the central section, z G [—8.0, 8.0] of the full field, as described in the main text. (Right) Maximum 
current |J| (solid line) and total kinetic energy (dashed line) and maximum Lorentz force (dot-dashed line) 
over time for the same field. 

• The technique used to transfer the initial state between codes has perturbed the magnetic field 
further from force- free. 

At the beginning of Sec. [3] we gave one suggestion, that the formation of the current layers may be 
due to an ideal instability as evidenced by the lack of dependence of kinetic energy on resistivity (see 
Fig.jsJ. Here we seek to determine additional information that may identify the cause of the dynamical 
evolution. We have noted that the integrated parallel current is a good predictor of the location of 
the two initial current layers. The integrated parallel current is a global quantity and so the global 
structure of the field may play a key role. The global structure is also important in, for example, the 
kink instability where a magnetic field with a set number of turns per unit length becomes unstable as 
more turns are added by increasing the length of the system. Whilst the kink-instability as it is usually 
considered applies to a tube with a well-defined single axis a similar kink-like instability, dependent on 
the total twist of the system, may apply to our braided field. 

With these considerations in mind we examine the evolution of only the middle section of the initial 
state of the braided magnetic field that is, we cut-out the section of the field in the above described 
experiment that lies in ^ G [—8, 8], x, 2/ G [—6, 6] at t = 0. This field is inserted as an initial condition in 
a new run, now keeping the fiux fixed on z = ±8, the new upper and lower boundaries of the domain. To 
maintain consistency we use the same resolution in the horizontal direction 320^ and a similar effective 
resolution in the vertical direction, 128 cells over z G [—8,8]. 

In the evolution of this new 'middle-third' field we find the system adjusts from its initial condition 
(with zero plasma velocity) to an approximate equilibrium in which the current structure is qualitatively 
similar to that of the initial state. To illustrate, contours of current in the central plane (z = 0) are shown 
at t = 20 in Fig. M (right) which may be compared with Fig. p\ (upper left) where the corresponding 
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contours in the initial state, t = 0, are shown. The maximum current in the domain is shown as a sohd 
hne in the left-hand image of Fig. |9] (left-hand axis), the kinetic energy integrated over x,y ^ [—3,3], 
z G [— 8, 8] as a dashed line (left-hand axis) and the maximum Lorentz force within the domain as a dot- 
dashed line (right-hand axis). As in the evolution of the full field, the Lorentz force in this initial state is 
high as a result of the transfer between grids and decreases rapidly over the first few time units. However, 
the Lorentz force subsequently stabilises at a low value as the system readjusts to equilibrium. The result 
that the middle section of the braid alone is in an equilibrium state suggests that an instability is present 
in the full braided field (rather than a lack of equilibrium due to numerical artefacts). Furthermore the 
instability is of a long-wavelength type with the full structure of the braided field being key. 

Returning to the evidence of Figure [s] during the time t G [8, 12] where the instability is clear in the 
current growth there is only a very slight dependence of the maximum velocity within the domain on 
resistivity and no increase in kinetic energy. This effect may be due to the confining nature of the strong 
background field external to the braided field structure which results in a deflection of the outflowing 
plasma around the boundary of the braided field (see Fig. |7|. Prior to t 8 no dependence of the flow 
on resistivity is seen suggesting an ideal dynamics where the system is adjusting to the distance from 
equilibrium. The current growth in t G [8, 12] does have a clear separation according to resistivity 77. 
However the increase is linear (rather than exponential) suggesting a dominant non-linear phase and the 
growth is slower for higher resistivity suggesting the non-linear phase is damped by resistive effects. 

5 Conclusions 

A previous paper (Wilmot-Smith et al. 2009a) considered the ideal relaxation of a braided magnetic field 
towards a force-free equilibrium. Here we have taken the final state of that relaxation and used it as an 
initial condition for a full resistive MHD simulation. 

The braided field is not in a stable equilibrium; two thin current sheets form after a short time (around 
a quarter of the time for an Alfven wave to cross the numerical box in the vertical direction) . The linear 
rate of increase of current density and the maximum strength of the current is found to increase with 
decreasing resistivity although the evolution of the total kinetic energy in the domain is independent 
of resolution. We conclude that the instability is possibly an ideal one although the details of how it 
occurs have not been determined. The wavelength of the instability was tested by considering the MHD 
evolution of the middle third section of the braid alone and this new field was found to be stable. Hence 
a long- wavelength dependence is implied. 

The initial configuration contains many regions of high squashing factor, Q, corresponding to quasi- 
separatrix layers. Plasma flow across QSLs is often thought likely to lead to current sheet formation. The 
two current sheets that form here do align with two of the highest regions of Q although the remaining 
regions of high Q do not correspond to any particular current features. The locations of the two current 
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layers are well predicted by the peak values of the integrated parallel current in the initial state. In the 
central plane this quantity shows two clear peaks and it is at these regions that the two current layers 
later form. 

These two current layers correspond to reconnect ion sites. Perpendicular to the layers the magnetic 
field has an elliptic structure, an admissible property of three dimensional reconnection that has only 
recently been found. The flow about the reconnection sites is of an asymmetric stagnation type, although 
large-scale rotational flows dominate the global structure; the effect of reconnection is not a strong 
acceleration of the flow but a more subtle untwisting process leading to the change in magnetic field 
topology. 

The longer-term evolution of the system will be considered in a future paper. 

Appendix 

Equation ^ can be derived from an evolution for the vector potential of a frozen-in magnetic field: 

^ + V (V- A) - V X V X A = 0. (8) 

This equation is equivalent to the Lie-derivative for a differential one-form, a = Aidx^, associated with 
the vector field A (Hornig, 1997). Hence (|8| can be written as 

^+L^a = 0, (9) 

and we can express a{t = 0) = F''a{t) (Abraham et al 1988), or more conveniently a{t) = (F~-^)*a(t — 
0), where F : X — > x maps the initial to the final coordinates and the star indicates the pull-back 
operation. This last equation written in components of the vector potential is ([t]). 
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